17th Computational Fluid Dynamics Conference, June 6-9, Toronto, Canada 


AIAA-2005-5229 


A Study of Multigrid Preconditioners 
Using Eigensystem Analysis 


Thomas W. Roberts* and R. C. Swanson'*' 

Mail Stop 128, NASA Langley Research Center, Hampton, VA 23681-2199 


The convergence properties of numerical schemes for partial differential equations are 
studied by examining the eigensystem of the discrete operator. This method of analysis 
is very general, and allows the effects of boundary conditions and grid nonuniformities 
to be examined directly. Algorithms for the Laplace equation and a two equation model 
hyperbolic system are examined. 


I. Introduction 

One of the critical needs in computational fluid dynamics is the development of rapidly converging 
algorithms. Present Reynolds averaged Navier-Stokes (RANS) schemes have exceedingly poor convergence 
rates on realistic problems. Understanding the source of slow convergence should lead to improvements in 
the algorithms. Although many sources of poor convergence are understood in a general, qualitative way 
(e.g., grid stretching and skewness), the exact nature of these difficulties in practical problems is often a 
mystery. 

The most common tool for analyzing the stability and convergence of numerical methods is Fourier 
analysis, especially von Neumann stability analysis . 1 This is applicable to uniform, Cartesian grids and 
constant coefficient equations, and neglects the effects of boundary conditions. Despite these seemingly 
severe restrictions, Fourier analysis is remarkably useful and leads to predictions of behavior that frequently 
are good quantitatively as well as qualitatively, even for real, nonlinear problems. However, boundary 
conditions can have a profound effect on the performance of a numerical scheme, and such effects cannot be 
determined with Fourier analysis. Understanding more general difficulties, such as the effects of stretched and 
unstructured grids, grid cell aspect ratio, and variable coefficient equations, requires more general analysis 
tools. 

One approach is to look at the eigensystem of a flow solver. Eriksson and Rizzi 2 used this approach 
to examine the behavior of their multistage, central difference Euler solver. Linearizing about a partially 
converged state, they used Arnoldi’s method to extract the slowest converging eigenmodes. In this way they 
were able to examine in detail the persistent modes of the scheme. At the time of their work, computer 
speeds and memory sizes were such that it was extremely expensive to extract a handful of modes. Present 
day hardware, as well as improved algorithms for extracting the eigenmodes of large systems , 3 makes their 
approach more viable. Furthermore, if a Krylov method such as GMRES is used as the iterative method, 
finding the eigenvalues and eigenvectors is very little additional work and provides a good diagnostic ca- 
pability . 4 On the other hand, this approach can be too general: it is not amenable to studying numerical 
difficulties in isolation and to evaluating the effectiveness of different components of an algorithm. One 
would like to be able to study simpler model problems, which may be clumsy to handle with a full-blown 
flow solver. 

Current computer speed and memory, and the availability of software packages such as Matlab®, makes 
it possible to compute the complete eigensystems of moderately sized discretizations of partial differential 
equations. If the iterative method is to be used as a preconditioner for a Krylov method, knowledge of the 
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spectrum can allow one to make good estimates for the convergence of the Krylov iteration. 5 In the present 
work, we explicitly calculate the eigensystems of general model equations. We restrict ourselves to linear 
systems, but allow for nonuniform grids and arbitrary boundary conditions. 

In the next section, a general description of iterative methods for solving linear systems is presented. 
Although this material is not new, it provides a convenient framework for developing the analysis tools. 
Multigrid preconditioning is described in some detail as a prelude to our examples. In this paper we illustrate 
the use of the eigensystem analysis by examining the Laplace equation and a two equation model hyperbolic 
system. The results for the Laplace equation show significant differences with the usual Fourier analysis. 
Analysis is presented for an implicit scheme and a multistage scheme applied to the two equation system. 

II. Structure of iterative schemes 

All linear, two level iterative methods for the solution of the system Lq = f can be viewed as precondi- 
tioned iterative methods. This includes classical relaxation methods such as Jacobi and Gauss-Seidel, as well 
as explicit and implicit time marching methods. In this section, the connection between most of the common 
methods for solving systems of linear equations arising from discretizations of partial differential equations 
is made explicit. Only linear systems of equations are considered. Although we are interested in discrete 
approximations to partial differential equations, the discussion is quite general. What follows is not new, 
but the authors have not seen this idea presented in this fashion previously. Lomax and Steger 6 discussed 
iterative methods in computational fluid dynamics in a very general way, which comes close to what follows 
herein. However, they did not make explicit the connection between the power series approximation to the 
inverse of a linear operator and most of the common methods for solving the system of equations. 

A. Power series approximation to the inverse 

Consider the function y = 1 — x, where x is a real or complex scalar. When |cc| < 1, then the binomial 
theorem gives the power series representation of the inverse of y , 

-1 1 9 

y = = 1 + x + x H 

1 — x 

OO 

= X>‘- 

fc =0 

This notion can be extended to a linear operator in an m dimensional space. First, consider the m x m 
matrix A. If det(I — A) ^ 0, where I is the identity matrix, and ||A|| < 1, then 

( I - A) -1 = 1 + A + A 2 H 

OO 

= E Ak a) 

k—0 

This power series representation of the inverse of a linear operator is especially attractive when m is large 
and A is sparse. A direct inversion of / — A may not be feasible either in storage or time. Although A is 
sparse, the inverse (I — A)^ 1 generally is not. If the series in Eq. (1) converges rapidly, then the inverse can 
be computed in 0(m) operations. 

Consider the system Lq = f. To invert this system iteratively, one almost invariably is required to 
precondition the system. Let P be the preconditioner, which is an m x m, nonsingular matrix. The 
preconditioned system to be solved is PLq = Pf. By choosing the preconditioner such that ||7 — PL || < 1, 
then PL = I— (I— PL) has the form in Eq. (1) with A = I— PL. This leads to the power series representation 
ofL" 1 

OO 

L- 1 =^2{I - PL) k P. (2) 

fc= o 

The preconditioner P may be thought of as an approximate inverse to L. The closer P is to L , the more 
rapid will be the convergence of the iteration. The goal is to find a P that is not only a good approximation 
to A -1 , but is sparse and easy to evaluate. In practice the preconditioner P is not usually formed explicitly. 
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Let L n x = — PL) k P, the n-term approximation to the inverse of L. Successive approximations 

to L satisfy the recurrence equation 


L-'=P+{I-PL)L~\ (3) 

In practice, we are interested in the solution to the system Lq = f rather than L -1 itself. It is usual to 
compute the corresponding sequence of approximations to <7, 

Qn = L~ l f 

= Pf+(I-PL)L~ 1 _ 1 f 

= Pf + (I-PL)q n _ 1 . (4) 

It is seen that the solution q is a fixed point of the iteration (4), and that the sequence converges to the 
solution of the original problem. 

If a sequence of v preconditioners Pk is used to construct L” 1 , 

^n-u+k = P k + ~ PkL) L n u+k l , k = 1 . . . V, 

then we can introduce a preconditioner P v defined by the recurrence equation 

Pi = Pi, Pk = Pk + {I — PkL) Pk- 1, k = 2...u (5) 

which allows us to write 

l? = p v + (i-p v l) l-\. 

The sequence of v preconditioning stages can be interpreted as subiterations, or an inner cycle, and P v is 
the preconditioner of the outer iteration or cycle. This recursive structure can be used to build a complex 
preconditioner from a general sequence of simpler ones. For example, a multigrid preconditioner is built up 
from a sequence of fine grid smoothing and coarse grid correction preconditioning steps. 

Although we have treated only left preconditioning, the above discussion applies as well to right precon- 
ditioning, i.e. , LPP~ 1 q = /, where LP rather than PL is the preconditioned operator. 

B. Illustrative examples of common iterative methods 

Common relaxation schemes are constructed by splitting the operator L into a sum of the form L = M — N , 
where M is nonsingular. Then an iterative scheme is 

Mq n = f + Nq n _i. 

Familiar examples are to take M to be the diagonal elements of L, which yields Jacobi iteration, or to take M 
to be the lower triangular part of L, which yields lexicographic Gauss-Seidel. Solving this eciuation gives 

q n = M^ 1 f + M~ 1 Nq n -i 

= M- 1 f + M~ 1 (M-L)q n _ 1 

= M- 1 f+(l-M~ 1 L)q n _ 1 . (6) 

Comparing Eq. (6) to Eq. (4), it is immediately apparent that this is a preconditioned iterative scheme 
where M~ l is the preconditioner. 

If L is the steady state part of a time dependent operator, dtq + Lq = /, then the solution to Lq = f can 
be obtained by marching in time until the asymptotic state is reached. This is the most common solution 
method in computational fluid dynamics. If explicit, forward time integration is used, and q n is the solution 
at time t = nr, then the iteration is 

q n = rf + (I -rL)q n _ 1 . (7) 

This has the form of Eq. (4) with a scalar preconditioner r. If local time steps are used to accelerate 
convergence, then r is a diagonal matrix with the local time steps along the diagonal. Similarly, explicit 
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multistage schemes and two level implicit schemes can be put in the form of Eq. (4). For example, the three 
stage explicit scheme 

q {0) = q n - 1, 

q {1) = g(°) + a X T (/ - Lg(°)) , 

<? (2) = g (0) + a 2 T (/ - , 

<Z ( 3 W 0) + r(f-Lq^), 

q n =q (3) 

is equivalent to the preconditioned scheme (4) with P = tI — o 2 t 2 L + aia 2 T 3 L 2 . Finally, consider the 
implicit approximate factorization scheme 

(/ + tLi) (/ + tL 2 ) ( q n - q n - 1) = t (/ - Fg„_i) , 

where P = Pi + P 2 . The preconditioner of this system is P = (I + rP 2 ) 1 (J + tPi) 1 r. 


C. Multigrid preconditioning 

The preconditioners in subsection II. B. are straightforward to write down explicitly. In more general cases 
the recursive expression Eq. (5) may be used to build an expression for the preconditioner corresponding to 
any arbitrary scheme in a straightforward way. We use this to construct the preconditioner corresponding 
to a general multigrid operator. 

A multigrid iteration consists of two parts, fine grid smoothing and a coarse grid correction. It is typical 
to parameterize a multigrid cycle by three numbers: iq, the number of smoothing iterations before the coarse 
grid correction; ^ 2 , the number of smoothing iterations after the coarse grid correction; and the cycle index 7, 
which is the number of coarse grid iterations in the cycle. A cycle index of one corresponds to a E-cycle, two 
corresponds to a IT-cycle. 

We wish to construct the preconditioner P m of a general multigrid cycle. The pre- and post-smoothing 
iterations have the same recursive structure shown in Eq. (4). Let P s be the preconditioner corresponding 
to the the particular smoother being used. We start to build P m by the iteration 

P m =0 

for k = 1 to q do P m <— P s + (I — P S L) P m . (8) 

Next a coarse grid correction is determined. Let P c be the preconditioner corresponding to a coarse grid 
update, i.e., the approximate inverse to the coarse grid operator L c . Let 1Z be the restriction operator 
(transfers data from a fine grid to a coarse grid) and V be the prolongation operator (transfers data from a 
coarse grid to a fine grid). The contribution to the preconditioner, P 7 , for a cycle index 7 is 

Pc = 0 

for k = 1 to 7 do P c <— P c + (/ — P C L C ) P c . 
p 7 = vP c n. 

This is used to update P m : 

Pm *■ — P'y + (I ~ P~/P) P m - 

After the coarse grid correction is computed, jq smoothing sweeps finish the cycle: 

for k = 1 to V 2 do P m <— P s + (I — P S L) P m . 

One multigrid cycle may be written q n = P m f + (I — P m L ) q n -\. 

The preconditioner P m can be written explicitly by expanding the terms in Eqs 


(9) 

(10) 

( 11 ) 

.(8) to (11), which yields 


V 2 — 1 

Pm = £ (/ - PsL) k P S + (I~ P s Lf 

k—0 


fv i — l 


f 7— 1 


v-i-1 


Y, (I ~ P s L) k P s + v (Y(I - PcLc) k Pc) n ( I - ( I~P s L) k P s 


( 12 ) 


, k—0 


\k = 0 


k—0 
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Because of the recursive nature of the multigrid cycle, the computation of the coarse grid preconditioner P c 
is done in the same way as P m . 

The convergence of the multigrid scheme can be analyzed by determining the spectrum of the amplification 
matrix G m = I — P m L. The asymptotic convergence rate is the modulus of the largest eigenvalue of G m . 


III. Model equation analysis 

The algorithm for the computation of P m and G m has been programmed in Matlab® for two model 
equations: the Laplace equation and a two equation hyperbolic system. The discrete operator L and the 
preconditioners P are represented explicitly as matrices. The spectrum of G m is computed using the built-in 
function eig. 


A. Two dimensional Laplace equation 

The first example we consider is the Laplace equation in two space dimensions, 


d x q + dyQ = 0 


(13) 


on the unit square with Dirichlet boundary conditions. Equation (13) is discretized with the standard five 
point central difference stencil. The difference equation is solved using multigrid with lexicographic point 
Gauss-Seidel as the smoother. As a point of reference, we initially take a uniform discretization in the 
x and y directions, so that a comparison with Fourier analysis may be made. 

The convergence of this scheme is usually estimated by smoothing analysis. 7 Letting 0 X and 9 y be 
the Fourier modes in the x and y directions, respectively, the amplification factor g of one sweep of the 
Gauss-Seidel smoother is 


g(0x,0 y ) = 1 - 


2 (cos 9 X — 1) + 2 (cos 9 y 

g — |— g 


1) 


(14) 


The modulus of the amplification factor \g\ <1/2 in the range n/2 < 9 X < 37r/2, n/2 < 9 y < 37t/2. 
This is the smoothing rate. Because the coarse grid is effective on the long wavelength error, the multigrid 
asymptotic convergence rate is estimated as the smoothing rate to the power v\ + ^ 2 , the total number of 
fine grid relaxation sweeps per cycle. 

A comparison of the spectra for lexicographic Gauss-Seidel relaxation and multigrid with the result of 
Fourier analysis is shown in Figure 1. Equation (13) is discretized on a uniform, 33 x 33 grid. The spectrum 
of the matrix operator for Gauss-Seidel with periodic boundary conditions is seen to compare extremely well 
with the result of Fourier analysis, as it should. When Dirichlet boundary conditions are introduced, however, 
the spectrum of the operator changes quite dramatically. The eigenvalues are more tightly clustered about 
the origin, and the slowly converging modes have collapsed to the positive real axis. This is qualitatively 
much different than the picture given by the Fourier analysis. When Gauss-Seidel is used as a smoother for 
a V(1,0) multigrid cycle (W (1,0) means v \ = 1, za> = 0) on five grids, all the eigenvalues are now well inside 
the limit given by the smoothing analysis. Another effect of the coarse grid correction is that the eigenvalues 
are not quite so tightly clustered near the origin. 

To confirm the validity of the analysis, Table 1 shows the convergence rates predicted from consideration 
of the maximum eigenvalue, and the actual convergence rates of the iterative solver for homogeneous Dirichlet 
boundary conditions and random initial conditions. The rates agree extremely well. We have also considered 
a four stage multistage scheme as the smoother, treating the Laplace operator as the steady state asymptote 
of the heat equation. A comparison of the convergence rates for this scheme with the predicted values is 
shown in Table 2. The standard multistage coefficients were chosen, and the time step was at the stability 
limit given by von Neumann analysis. Again, the agreement between the actual and predicted convergence 
rates is excellent. 

One of the advantages of this analysis is that we can examine the effect of variable coefficients, such 
as arise from grid stretching. Equation (13) is discretized on a unit square with uniform spacing in the 
^-direction and a stretching rate of 1.2, a severe stretching, in the y-direction. The aspect ratio of the grid 
cells, h y /h x , varies from 0.0188 at the lower boundary to 5.349 at the upper boundary. This is a case that 
cannot be treated by Fourier analysis. The effect of grid stretching is shown in Figure 2, where both the 
Gauss-Seidel smoother alone and a V(1,0) cycle are compared to the uniform grid case. The degradation of 
the convergence rate is apparent, but about 90% of the eigenvalues have a modulus less than 1/2. Except for 
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Table 1. Predicted and actual multigrid convergence rates for the solution of the Laplace equation on a uni- 
form 33 x 33 grid using lexicographic Gauss-Seidel smoothing and homogeneous Dirichlet boundary conditions 


Grids 

V(l, 0) cycle 
Actual Predicted 

W (1,0) cycle 
Actual Predicted 

1 

0.9901 

0.9904 

- 

- 

2 

0.9530 

0.9530 

0.9170 

0.9170 

3 

0.8191 

0.8191 

0.5006 

0.5006 

4 

0.4658 

0.4658 

0.3033 

0.3016 

5 

0.3314 

0.3318 

0.2975 

0.3016 


Table 2. Predicted and actual multigrid convergence rates for the solution of the Laplace equation on a 
uniform 33 x 33 grid using four stage multistage smoothing and homogeneous Dirichlet boundary conditions 


Grids 

V(l, 0) cycle 
Actual Predicted 

TV (1,0) cycle 
Actual Predicted 

1 

0.9952 

0.9952 

- 

- 

2 

0.9764 

0.9764 

0.9579 

0.9579 

3 

0.9074 

0.9074 

0.7153 

0.7153 

4 

0.7116 

0.7116 

0.6038 

0.6053 

5 

0.6028 

0.6058 

0.6033 

0.6053 


the persistent modes lying along the positive real axis, the spectra are qualitatively similar to their uniform 
grid counterparts. 

B. Two equation hyperbolic system 

The second example considered is the model two equation system 


dtu — nd x u — d y v = 0, (15a) 

dtv + d x v — d y u = 0. (15b) 

This is a particularly simple hyperbolic model system. When the parameter k = 1, steady solutions to the 
system (15) satisfy the Cauchy- Riemann equations. If k is taken to be the Prandtl-Glauert factor, steady 
solutions to Eq. (15) can be interpreted as the velocity components of a small perturbation compressible 
flow. 8. 

To solve this system, we initially consider the implicit scheme described by Roberts and Warren 8,9 to 
discretize the system. First order upwind differencing is used on the left hand side, and second order upwind 
biased differencing is used on the right hand side. Using Fourier analysis, Roberts and Warren showed that 
if a direct inversion of the left hand side operator is performed, all wavelengths (except for 6 X = 0, 6 y = 0) 
have an amplification factor of less than 1/2. In practice, the left hand side is not directly inverted, but must 
be solved approximately. Red-black Gauss-Seidel iteration is used to approximately invert the first order 
left hand side operator. This is taken as the smoother for a multigrid cycle for the full system. The number 
of Gauss-Seidel subiterations of the first order system is an additional parameter of the iterative scheme. 

Figure 3 shows a comparison of spectra for Eq. (15) with n = 1 on a uniform 32 x 32 grid. Ten red- 
black subiterations were used to invert the first order operator and the CFL number was 10 6 * . The upper 
left plot shows the eigenvalues for a direct inversion of the left hand side. The dashed circle shows the 
maximum amplification factor given by Roberts and Warren. 8 On the upper right, the spectrum for ten 
red-black iterations and no multigrid is shown. The shape is very close to that for a direct inversion, with 
the exception of a few outliers. Incorporating this smoother into a five level V(1,0) multigrid cycle gives 

a However, Eq. (15) is not the correct system of unsteady small perturbation compressible flow. 
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Table 3. Predicted and actual multigrid convergence rates for the solution of two equation hyperbolic system 
on a uniform 32 X 32 grid using ten red-black relaxation sweeps to invert the left hand side 


Grids 

V(1,0) cycle 
Actual Predicted 

W(l, 0) cycle 
Actual Predicted 

W(l, 1) cycle 
Actual Predicted 

1 

0.7210 

0.7201 

- 

- 

- 

- 

2 

0.6361 

0.6361 

0.6568 

0.6570 

0.2603 

0.2598 

3 

0.6683 

0.6689 

0.6660 

0.6658 

0.2595 

0.2598 

4 

0.7373 

0.7373 

0.6676 

0.6671 

0.2606 

0.2598 

5 

0.7295 

0.7295 

0.6657 

0.6663 

0.262 

0.2598 


the plot on the lower left. Although there is more clustering of the eigenvalues near the origin, a large 
number of them he outside the dashed circle, leading to no improvement in the asymptotic convergence rate. 
Finally, a direct solve of the coarse grid problem gives the two grid spectrum on the lower right, which differs 
qualitatively from the V-cycle spectrum. The two grid convergence rate is 0.5649, which is still not as good 
as single grid results with a direct inversion of the left hand side. 

A comparison of the convergence rates for the system using various numbers of grids and cycles is shown 
in Table 3. All results are for a uniform 32 x 32 grid, k = 1, a CFL number of 10 6 and ten red-black 
subiterations of the first order operator as the smoother. The observed convergence rates for an iterative 
solution of Eq. (15) with homogeneous characteristic boundary conditions and random initial conditions 
are shown for comparison. The comparison of the predicted and actual rates are in excellent agreement. 
Interestingly, the addition of more grid levels does not improve the convergence rate. In fact, for the V(1,0) 
cycle the convergence deteriorates when more than one coarse grid level is used. The W(l,l) cycle is 
completely insensitive to additional coarse grid levels. 

We now consider a scheme recently introduced by Rossow. 10 This scheme uses an explicit Runge-Kutta 
method as a framework. An implicit stage is included in each Runge-Kutta stage in order to extend stability 
and appropriately treat geometrically stiff discrete systems, which generally occur in numerically resolving 
flows with boundary layers. The implicit operator corresponds to a first-order approximation of the governing 
flow equations, and the full operator is approximately inverted iteratively with Gauss-Seidel relaxation. In the 
present work, unless otherwise indicated, an upwind biased discretization of the explicit residual function was 
implemented as described by Roberts and Warren. For the analysis we used a 5-stage Runge-Kutta scheme 
and set the CFL number to 100. The coefficients for the Runge-Kutta scheme were taken from Van Leer et 
ah, 11 and they are as follows: 0.0695, 0.1602, 0.2898, 0.5060, 1. Lexicographic (Lex), Red-Black (R-B), and 
symmetric (Sym) Gauss-Seidel (G-S) iterations were considered for the implicit operator inversion. 

In the the upper plots of Figure 4 spectra for the Runge-Kutta scheme applied to the hyperbolic system 
of Eq. (15) on the uniform 32 x 32 grid are shown. For comparison purposes consider also the corresponding 
spectra for the implicit scheme in Figure 3. Both schemes use red-black Gauss-Seidel iterations for ap- 
proximately inverting the implicit operator. Here the implicit stage of the Runge-Kutta scheme uses three 
iterations. With the Runge-Kutta scheme as the smoother for the V(1,0) multigrid cycle, there is signifi- 
cantly stronger eigenvalue clustering, as evident in the upper right plot. This improved clustering is reflected 
in the convergence rates presented in Table 4. As exhibited with the implicit scheme, additional grid levels 
in the multigrid provide little or no improvement in the asymptotic convergence rate. 

The effect of mesh refinement on the uniform grid is also included in Figure 4. For the lower plots the 
analysis was applied on a 64 x 48 grid. The spectral patterns for the two grids are quite similar, differing 
primarily in density. There is the expected slowdown in single grid convergence rate (from 0.8002 to 0.8965). 
With the V(l, 0) multigrid cycle the convergence rate is 0.4508 on both meshes. 

Figure 5 shows the effect of grid aspect ratio alone on the spectra of the Runge-Kutta scheme with 
V(1,0) multigrid. In addition, the impact of different types of Gauss-Seidel relaxation is indicated. The 
number inside the parentheses denotes the number of iterations performed in each implicit stage. The 
aspect ratio is 100 throughout a 32 x 32 grid. With six iterations for approximate inversion of the implicit 
operator, the asymptotic convergence rates with lexicographic and red-black Gauss-Seidel are 0.6237 and 
0.6450, respectively. In the lower right plot we see that there is a strong effect on the spectrum clustering 
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Table 4. Predicted and actual multigrid convergence rates for the solution of two equation hyperbolic system 
on a uniform 32 X 32 grid using a five stage R-K scheme, with an implicit stage corresponding to each explicit 
stage. Three red-black relaxation sweeps are used to invert the implicit operator. 


Grids 

K(1,0) cycle 
Actual Predicted 

W(l, 0) cycle 
Actual Predicted 

W{ 1, 1) cycle 
Actual Predicted 

1 

0.8000 

0.8002 

- 

- 

- 

- 

2 

0.5237 

0.5240 

0.4522 

0.4522 

0.2039 

0.2038 

3 

0.4507 

0.4508 

0.4524 

0.4524 

0.2038 

0.2039 

4 

0.4504 

0.4508 

0.4524 

0.4524 

0.2040 

0.2039 

5 

0.4508 

0.4508 

0.4524 

0.4524 

0.2039 

0.2039 


with symmetric Gauss-Seidel, resulting in a convergence rate of 0.4579. Although not shown, even an aspect 
ratio of 1000 resulted in the same rate of convergence. This suggests that there is at most a weak sensitivity 
to aspect ratio. 

In general a good preconditioner needs to be effective not only on a grid with aspect ratio but also a grid 
with stretching. While Fourier analysis is not appropriate for stretched grids, the matrix analysis discussed 
herein can easily handle general stretchings in all coordinate directions. We now consider the effect of grid 
stretching on the performance of the Runge-Kutta scheme. For this purpose a 32 x 32 affine grid with 
variable spacing in the x- and y-directions was generated. The grid was defined on a rectangular domain 
with — 5 < x < 5 and 0 < y < 5. Points in the x-direction were distributed with a tangent stretching function 
on the intervals —5 < x < —1/2, 1/2 < x < —5, and a cosine spacing on the interval —1/2 < x < 1/2. In 
the y-direction the grid was stretched with a rate of 1.2 on the interval 0 < y < 5. The grid aspect ratio, 
h x /h y , varies from about 700 to about 1/90 over the computational domain. 

The spectra for the Runge-Kutta scheme applied to the stretched grid are displayed in Figure 6. The 
implicit operator was either approximately inverted with red-black Gauss-Seidel or symmetric Gauss-Seidel. 
For both types of Gauss-Seidel the same computational work was done. There is stronger eigenvalue cluster- 
ing and fewer outliers with the symmetric iteration. Nevertheless, with the K(1,0) multigrid cycle there is 
good clustering for both red-black and symmetric iteration. This suggests that with multigrid both versions 
of the Runge-Kutta scheme would be amenable to a Krylov method (e.g., GMRES). As smoothers for the 
multigrid, the two versions of the Runge-Kutta scheme have nearly the same maximum eigenvalue in modu- 
lus (0.7259 for red-black iteration and 0.7594 for symmetric iteration). To verify these predicted convergence 
rates, calculations were performed on the same grid used in the analysis for inviscid flow over a parabolic 
bump. A perturbation boundary condition was imposed at the solid boundary and the initial solution was 
set to zero. The convergence rates for the computations with the two variations of the Runge-Kutta scheme 
were essentially the same as the predicted rates. 

So far we have examined the performance of the Runge-Kutta scheme with only the wide discretization 
stencil (21 points in two dimensions) of reference. 8 In the application of the Runge-Kutta scheme Rossow 
used a narrow second order upwind differencing stencil (9 points). A question arises concerning the effect 
of the difference stencil type, especially since Roberts and Warren showed that a first-order implicit scheme 
with a narrow pure upwind second order stencil for the right-hand side has poor damping properties. In 
Figure 7 the spectra of the Runge-Kutta scheme for the aspect ratio(AR) of 100 grid and the stretched grid 
are displayed. Symmetric Gauss-Seidel was used to invert the implicit operator in all cases. The upper plots 
in the figure are for the AR = 100 grid, and the lower plots are for the stretched grid. In comparing the 
upper right plot with the lower right plot of Figure 5, one can see that the eigenvalue patterns for the two 
stencils are different. Nevertheless, for each stencil all the eigenvalues lie within the radius =1/2 circle. 
There is somewhat faster K(l, 0) multigrid convergence with the narrow stencil (rate = 0.3876) than with 
the wide stencil (rate = 0.4579). For the stretched grid the differences between the eigenvalue patterns for 
the narrow and wide stencils can be seen by comparing the lower plots of Figures 6 and 7. The convergence 
rates corresponding to the narrow stencil are again somewhat faster. With the single grid the narrow stencil 
rate is 0.8538, and the wide stencil rate is 0.8735. For the K(1,0) multigrid the convergence rates with the 
narrow and wide stencils are 0.6900 and 0.7594, respectively. 
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IV. Summary 


It has been shown that two level iterative schemes for the solution of a linear system Lq = f can be 
interpreted as solution to the preconditioned system PLq = P f by a power series expansion for (PL) -1 . This 
includes not only classical relaxation methods, but also two level time marching schemes for the unsteady 
equation dtq + Lq = f. Using this common structure leads to a simple algorithm for the computation of the 
preconditioner and amplification matrix. The spectrum of the latter can be efficiently found for moderately 
sized problems. This provides a very general and flexible way to study the convergence properties of numerical 
methods for partial differential equations. 

The two model equations presented here are illustrative of the use of the method. The authors are 
particularly interested in using the eigensystem analysis to address the convergence difficulties of the Navier- 
Stokes equations on highly stretched grids. The initial steps in extending the analysis to the compressible 
Navier-Stokes equations have been taken. Once this extension is completed we plan to consider the analysis 
for unstructured grids, to which Fourier methods do not apply. 
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Gauss-Seidel, Fourier analysis 



Gauss-Seidel, periodic boundary 



Gauss-Seidel, Dirichlet boundary 



V(1 ,0), Dirichlet boundary 



Figure 1. Spectra of preconditioners for the Laplace equation on a 33 x 33 grid. Solid line is the unit circle, 
dashed line is the smoothing rate given by Fourier analysis. 
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Gauss-Seidel, uniform grid 


V(1 ,0), uniform grid 




Gauss-Seidel, stretched grid 


V(1 ,0), stretched grid 




Figure 2. Spectra of preconditioners for the Laplace equation on a uniform and a stretched 33 x 33 grid. Solid 
line is the unit circle, dashed line is the smoothing rate given by Fourier analysis. 
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V(1,0), 10 subiterations 


Two-grid (1,0), 10 subiterations 




Figure 3. Spectra for the two equation system on a uniform 32 x 32 grid. Solid line is the unit circle, dashed 
line is amplification factor for a direct inversion of the first order left hand side operator. 
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Runge-Kutta (R-B G-S, 3) 


Runge-Kutta (R-B G-S, 3), V(1 ,0) 
32x32 32x32 



Runge-Kutta (R-B G-S, 3) 

64x48 



Runge-Kutta (R-B G-S, 3), V(1 ,0) 
64 x 48 



Figure 4. Spectra for the two equation system on uniform 32 x 32 and 64 x 48 grids. Solid line is the unit circle, 
dashed line corresponds to amplification factor of 0.5. 
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Runge-Kutta (Lex G-S, 6), V(1 ,0) 



Runge-Kutta (Sym G-S, 3), V(1 ,0) 



Figure 5 . Spectra for the two equation system on a uniform, aspect ratio of 100 , 32 x 32 grid. All cases for 
1/(1, 0) multigrid cycle. Solid line is the unit circle, dashed line corresponds to amplification factor of 0 . 5 . 
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Runge-Kutta (R-B G-S, 6) 



Runge-Kutta (R-B G-S, 6), V(1 ,0) 



Runge-Kutta (Sym G-S, 3) 



Runge-Kutta (Sym G-S, 3), V(1 ,0) 



Figure 6. Spectra for the two equation system on a stretched 32x32 grid. Grid stretched in x- and y-directions. 
Solid line is the unit circle, dashed line corresponds to amplification factor of 0.5. 
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Runge-Kutta (Sym G-S, 3), V(1 ,0) 


Runge-Kutta (Sym G-S, 3) 
AR = 100 



Runge-Kutta (Sym G-S, 3) Runge-Kutta (Sym G-S, 3), V(1 ,0) 

Stretched 





Figure 7. Spectra for the two equation system on a 32 x 32 grid. Standard nine point upwind differencing used 
for discretization. Top two plots for aspect ratio AR = 100 grid; bottom two plots for stretched grid. Solid 
line is the unit circle, dashed line corresponds to amplification factor of 0.5. 
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